library(survival)
library(stargazer)
library(mlogit)
library(sandwich)
library(lmtest)
library(sandwich)
library(lmtest)

clrobustse <- function(fit.model, clusterid) {
  rank=fit.model$rank
  N.obs <- length(clusterid)            
  N.clust <- length(unique(clusterid))  
  dfc <- N.clust/(N.clust-1)                    
  vcv <- vcov(fit.model)
  estfn <- estfun(fit.model)
  uj <- apply(estfn, 2, function(x) tapply(x, clusterid, sum))
  N.VCV <- N.obs * vcv
  ujuj.nobs  <- crossprod(uj)/N.obs
  vcovCL <- dfc*(1/N.obs * (N.VCV %*% ujuj.nobs %*% N.VCV))
  coeftest(fit.model, vcov=vcovCL)
}

cl.mlogit <- function(fm, cluster){
  
  M <- length(unique(cluster))
  N <- length(cluster)
  K <- length(coefficients(fm))
  # edit 6/22/2015: change dfc
  dfc <- (M/(M-1))*((N-1)/(N-K))
  # dfc <- (M/(M-1))
  uj <- apply(estfun(fm),2, function(x) tapply(x, cluster, sum));
  vcovCL <- dfc*sandwich(fm, meat.=crossprod(uj)/N)
  coeftest(fm, vcovCL) }

######################################
#### Reanalysis of Kreutz (2010) #####
######################################

## Loading replication data
kreutz <- read.csv("Kreutz_reanalysis.csv")


### Exact Replication of Kreutz 2010 (Table III)

K1_repl <- glm(formula = recur ~ vic + peace + part + ethnic + totalg + lnbds + lndur + infant_lag + demo_lag + elf +
             pcedur_yr + year + X_spline1 + X_spline2 + X_spline3, data = kreutz, family = "binomial")
K1_est_data <- na.omit(kreutz[,c("country", "case_id", names(K1_repl$coefficients)[-1])])
K1_est_test <- clrobustse(K1_repl, as.character(K1_est_data$country))

K2_repl <- glm(formula = recur ~ govwin + rebwin + peace + cease + pko + eth_rev + sec_ + lnbds + lndur +
                 perc_mil + infant_lag + gdppc_lag + demo_lag + lnpop + pcedur_yr, data = kreutz, family = "binomial")
K2_est_data <- na.omit(kreutz[,c("country", "case_id", names(K2_repl$coefficients)[-1])])
K2_est_test <- clrobustse(K2_repl, as.character(K2_est_data$country))

### Reanalysis based on repeated/recast conflict

## Guinea and Sudan are miscoded in Kreutz 2010. Excluding these cases.
dat1a <- subset(kreutz, kreutz$case_id != 197 & kreutz$case_id != 261)

casesplit <- split(dat1a, dat1a$case_id)

out <- data.frame()
for (i in 1:length(casesplit)){
  x <- casesplit[[i]]
  x <- x[order(x$year),]
  x$t_time <- 1:(nrow(x))
  out <- rbind(out, x)
}

## If both former and new rebels are fighting in the recurred war, it is coded as "repeated"
out$state1 <- out$recur * ifelse(out$status1 == 0, 0, ifelse(out$status1 == 3, 2, 1))
out$state2 <- out$recur * ifelse(out$status2 == 0, 0, ifelse(out$status2 == 3, 2, 1))

out$status <- out$state2
out$pid <- out$case_id
out$state <- ifelse(out$status2 == 0, 0, ifelse(out$status2 == 3, 2, 1))
out$event <- out$recur


# K_Reanalysis_model1 (Walter 2004)

K_reanal1_data <- na.omit(out[,c("status", "event", "vic", "peace", "part", "ethnic", "totalg", "lnbds", "lndur", "infant_lag", "demo_lag", "elf", "pcedur_yr", "year",
                      "X_spline1", "X_spline2", "X_spline3", "country", "case_id")])
K_reanal1_data$country <- as.character(K_reanal1_data$country)

ml_data1 <- mlogit.data(K_reanal1_data, shape = "wide", choice = "status")

K_res1 <- mlogit(formula = status ~ 0|vic + peace + part + ethnic + totalg + lnbds + lndur + infant_lag + demo_lag + elf +
                  pcedur_yr + year + X_spline1 + X_spline2 + X_spline3|0, data = ml_data1)

K_test1 <- cl.mlogit(K_res1, K_reanal1_data$country)

## Number of units and events
length(unique(as.character(K_reanal1_data$case_id)))
sum(K_reanal1_data$status == 1)
sum(K_reanal1_data$status == 2)

# model2 (Quinn et al. 2007)

K_reanal2_data <- na.omit(out[,c("status", "govwin", "rebwin", "peace", "cease", "pko", "eth_rev", "sec_", "lnbds", "lndur",
                                 "perc_mil", "infant_lag", "gdppc_lag", "demo_lag", "lnpop", "pcedur_yr", "country", "case_id")])
K_reanal2_data$country <- as.character(K_reanal2_data$country)

ml_data2 <- mlogit.data(K_reanal2_data, shape = "wide", choice = "status")

K_res2 <- mlogit(formula = status ~ 0|govwin + rebwin + peace + cease + pko + eth_rev + sec_ + lnbds + lndur +
         perc_mil + infant_lag + gdppc_lag + demo_lag + lnpop + pcedur_yr|0, data = ml_data2)
K_test2 <- cl.mlogit(K_res2, K_reanal2_data$country)
length(unique(as.character(K_reanal2_data$case_id)))
sum(K_reanal2_data$status == 1)
sum(K_reanal2_data$status == 2)

## Arranging the results for the table output
K_10 <- K1_est_test
K_11 <- K_test1[seq(1, nrow(K_test1), 2),]
row.names(K_11) <- rownames(K_10)
K_12 <- K_test1[seq(2, nrow(K_test1), 2),]
row.names(K_12) <- rownames(K_10)

K_20 <- K2_est_test
K_21 <- K_test2[seq(1, nrow(K_test2), 2),]
row.names(K_21) <- rownames(K_20)
K_22 <- K_test2[seq(2, nrow(K_test2), 2),]
row.names(K_22) <- rownames(K_20)

### Table 1

K_reanal_table <- stargazer(K_10, K_10, K_10, K_20, K_20, K_20, type = "latex",
          coef = list(K_10[,1], K_11[,1], K_12[,1], K_20[,1], K_21[,1], K_22[,1]),
          se = list(K_10[,2], K_11[,2], K_12[,2], K_20[,2], K_21[,2], K_22[,2]),
          title = "Results from Reanalysis of Kreutz (2010)", style = "aer",
          digits = 2,
          no.space = T,
          order = c(1,2,3,4,8,5,9,6,10,7,11,12,13,14,15,16,17,18,19,
                    20,21,22,23,24,25),
          column.labels = c("Replication", "Repeated", "Recast", "Replication",
                            "Repeated", "Recast"),
          covariate.labels = c("Victory", "Government victory", "Rebel victory", "Peace agreement",
                               "Ceasefire", "Partition", "Peacekeepers",
                               "Ethnic mobilization", "Ethnic revolution",
                               "Total goals", "Secessionist",
                               "Battle deaths (ln)", "Duration (ln)", "Army size (% of pop.)",
                               "Infant mortality (lag)", "GDP/capita (lag)", "Democracy (lag)",
                               "ELF", "Population (ln)", "Peace years", "Year", "Spline (1)",
                               "Spline (2)", "Spline(3)"))
## To save
#cat(K_reanal_table, sep = '\n', file = 'K_table.tex')

############################################################
#### Reanalysis of Doyle and Sambanis 2006 (Appendix B) ####
############################################################

ds <- read.csv("D&S2006_reanalysis.csv")

sv <- Surv(ds$duration, event = ifelse(ds$peacend == 1, 1, 0)) 

## Exact Replication
dscox1 <- coxph(sv ~ wartype + logcost + factnum + transpop +
                  idev1 + gdpgrofl + isxp2 + unintrvn + treaty +
                  cluster(clust2), ds, method = "breslow")


sv1 <- Surv(ds$duration, event = ifelse(ds$peacend2 == 1, 1, 0))
sv2 <- Surv(ds$duration, event = ifelse(ds$peacend2 == 2, 1, 0))

## Reanalysis based on repeated/recast conflict
dscox2 <- coxph(sv1 ~ wartype + logcost + factnum + transpop +
                  idev1 + gdpgrofl + isxp2 + unintrvn + treaty + 
                  cluster(clust2), ds, method = "breslow")
dscox3 <- coxph(sv2 ~ wartype + logcost + factnum + transpop +
                  idev1 + gdpgrofl + isxp2 + unintrvn + treaty +
                  cluster(clust2), ds, method = "breslow")

### Table 4
stargazer(dscox1, dscox2, dscox3,
          omit.stat = c("rsq", "max.rsq", "wald", "logrank"), style = "ajps")
